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Abstract 

We present necessary and sufficient optimality conditions for finite time optimal con¬ 
trol problems for a class of hybrid systems described by linear complementarity models. 
Although these optimal control problems are difficult in general due to the presence of 
complementarity constraints, we provide a set of structural assumptions ensuring that the 
tangent cone of the constraints possesses geometric regularity properties. These imply that 
the classical Karush-Kuhn-Tucker conditions of nonlinear programming theory are both 
necessary and sufficient for local optimality, which is not the case for general mathematical 
programs with complementarity constraints. We also present sufficient conditions for global 
optimality. 

We proceed to show that the dynamics of every continuous piecewise affine system can 
be written as the optimizer of a mathematical program which results in a linear complemen¬ 
tarity model satisfying our structural assumptions. Hence, our stationarity results apply to a 
large class of hybrid systems with piecewise affine dynamics. We present simulation results 
showing the substantial benefits possible from using a nonlinear programming approach to 
the optimal control problem with complementarity constraints instead of a more traditional 
mixed-integer formulation. 

Many dynamical system have both continuous and switching, i.e. hybrid, dynamics [26]. 
Common examples include mechanical systems with impact [6], traction models for cars [4, 
16], and biological systems [43]. The class of piecewise affine (PWA) systems is particularly 
important because it is a natural generalization of linear systems and hence conceptually easy 
to understand and use for modeling [52]. It has been shown in [30] that PWA systems are 
equivalent to a number of other hybrid system model classes such as linear complementarity 
(LC) models [29], mixed-logical dynamical (MLD) [1], and max-min-plus-scaling models [15] 
under certain technical assumptions. Some tools and results developed for one system class 
may, accordingly, be applied to another, equivalent system class where appropriate. Finite 
horizon optimal control problems for discrete time hybrid system models have for example 
been studied in [1] (PWA / MLD models) and [13, 14] (max-min-plus-scaling models). 
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Here we consider hybrid dynamical systems in LC form: 

x~^ = Ax + BuU + BwW + c, (la) 

0 < EyjW + ExX + EuU + e _L w > 0, (lb) 


where x G is the current system state and x~^ G the successor state resulting from 
applying a control input u G M”“. Since both sides of (lb) are constrained to be non-negative, 
the orthogonality relation _L requires that every component of the complementarity variable 
w G is constrained to be zero if the corresponding component on the left-hand side of (lb) 
is non-zero, and vice-versa. While the complementarity variables w can by ascribed physical 
meaning with respect to the modeled system in some cases, see [3, Ch. 11], [19], or [53], 
we will consider the general case where they are treated as arbitrary auxiliary variables. The 
other matrices in (1) are assumed to have compatible dimensions, in particular E^j G 
is a square matrix. 


For the hybrid system (1) we consider finite horizon optimal control problems: 

7V-1 

min iNixN) + ik{xk,Uk) 

11.X.W f ^ 


k=0 


s.t. Xk +1 = Axk + BuUk + BxuWk + C 

0 < E^ijWk + ExXk + EuUk + e _L Wk > 0 


(2a) 

(2b) 

(2c) 


Here, Xk G denotes the predicted system state k timesteps into the future and the 
constraints have to be satisfied for all A: = 0,..., N—1. We denote by u := (uq, ..., the 

collection of the control inputs to be chosen, by x := (xq, ... ,xn) the resulting state trajectory 
starting from the initial state xq = x, and by w := {wq, ... ,wn-i) the corresponding 
trajectory of complementarity variables. We assume the typical case of convex stage costs 
ik- 1^”“’ X —)■ M and terminal cost —)■ M that are independent of Wk- The 

constraints (2b) and (2c) have to hold for all k = 0,..., N — 1. The variables Wk do not 
appear in the the cost function (2a) because we consider them to be auxiliary variables without 
physical meaning. 

We deliberately omit constraints on the state and control input from the optimal control 
problem (2) for clarity of exposition. However, all of the results in this paper extend to linear 
state-input constraints over the horizon, i.e. {xk,Uk) G T^ for polytopes T^ C x M"^". 

Due to the complementarity constraints (2c) that are part of the LC dynamics, the opti¬ 
mal control problem (2) is called a mathematical program with complementarity constraints 
(MPCC) [42]. For an overview of MPCC solution methods and different applications see [3, 
46]. MPCCs are nonlinear and nonconvex optimization problems (i.e. the complementarity 
constraint (2c) can be modeled as a bilinear (in-)equality) and are generally difficult to solve. 
The main reason for this is that the classical Mangasarian-Fromovitz constraint qualification 
(MFCQ) of nonlinear programming (NLP) is violated at every feasible point [8]. MFCQ is 
typically assumed in NLP theory to enable the use of Karush-Kuhn-Tucker (KKT) conditions 
to characterize optimality. It is equivalent to the compactness of the set of Lagrange multipliers 
at a local optimum [25]. 
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In this paper we make a number of structural assumptions on (1) to focus on a particular 
subclass of LC models. Under these assumptions we can prove strong stationarity conditions 
for the optimal control MPCC (2). In particular, we will show that the classical KKT conditions 
from NLP theory are both necessary and sufficient for optimality in (2). This is generally not 
the case for MPCCs where KKT conditions cannot be expected to hold and specialized solution 
methods are often used [42]. 

We then show that our results for (2) cover optimal control problems for arbitrary continuous 
PWA systems. In particular, we exploit the fact that continuous PWA functions can be written 
as the difference of two convex PWA functions [33, 40] to represent PWA system dynamics 
as the solution to a convex parametric quadratic program (PQP). Further manipulations result 
in a special LC model (1) which satisfies the assumptions we require on (1) for our strong 
stationarity results. We also present sufficient conditions for a given local optimum to be the 
only local optimum within a neighborhood around it or even globally optimal for (2). 

In addition to an illustrative example we present numerical results that indicate solving (2) as 
an NLP can be computationally advantageous compared to more traditional mixed-integer 
approaches to solving optimal control problems for PWA systems [1]. Treating (2) as a 
general NLP instead of an MPCC is only possible due to the strong stationarity results in 
this paper since typical NLP algorithms attempt to find a solution to the KKT conditions 
which, as mentioned above, are not suitable optimality conditions for MPCCs. Solving (2) 
as an NLP leads to significantly shorter computation times and much better scaling in the 
problem dimensions than an MIP approach and often yields globally optimal solutions in 
numerical experiments. While no a-priori guarantees exist for this, simulation studies for 
randomly generated systems outline the possible benefits of this approach. Since optimal 
control problems for hybrid systems are NP hard in general no such guarantees should be 
expected [10]. 

Notation: We use superscripts in square brackets to indicate elements of vectors and vector 
valued functions and matrices, e.g. is the j* element of / and is the row of M. 
We use the symbols A and V to denote the logical operations “and” and “or”, respectively. 
The Kronecker product between two matrices A and B is denoted A®B. The notation x Ay 
indicates that the two vector x,y are orthogonal, i.e. x~^y = 0. A bold 1 denotes a 

vector (or matrix) of ones and I an identity matrix. In case the dimensions are not clear from 
context we indicate them using subscripts, e.g. Inxm £ is an n-by-m matrix of ones. 

A positive (semi-)definiteness condition on a symmetric matrix M = M~^ is denoted M A 0 
(M A 0). While H-H denotes an arbitrary norm on the respective space, we denote a weighted 
Euclidean norm using a symmetric matrix Q >- 0 as H^Hq := y/ x~^Qx. The sign of a scalar 
X is given by sgn(x). 

The nullspace of a matrix M is given by M{M). The interior of a set U C is denoted 
int(U) and the cardinality of a set 7 C N by |7|. We denote the set of real non-negative 
numbers M_|_ and the index set {l,...,n} by N^. All vector-valued inequalities are to be 
understood component-wise. 

The term mathematical program with equilibrium constraints (MPEG) is also used regularly 
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in the literature for MPCCs such as (2), although the two prohlem classes are not entirely 
equivalent [12]. We will use the term MPCC throughout the paper, except for some technical 
terms, e.g. MPEC-Ahadie constraint qualification, which have been firmly established in the 
literature using the MPEC-acronym. 

I. Standing assumptions and main result 

Throughout the paper we make a number of standing assumptions about the EC model (1). 
The first is a reasonable assumption to guarantee deterministic behavior: 

Assumption 1. Every complementarity variable w that solves {Vo) for a fixed {x,u) results 
in the same successor state x^. 

Eor an in-depth discussion of well-posedness of EC systems and related results we refer to [7]. 

The next is a technical assumption on the structure of the complementarity problem (lb): 

Assumption 2. There exist an nf, G N and scalars ,. .., G N that partition 
such that the complementarity problem (lb) can (possibly after a coordinate transformation) 
be decomposed into Uh independent complementarity problems 

Vi G Nnt: 0 < MiWi + CiX + DiU + 6* _L Wi > 0, (3) 

where Mi = MJ G is a symmetric positive semidefinite rank-one matrix, i.e. 

Vi G Nnt : Mi = rriimj A 0 

f?l 

for some rrii G This decomposition is minimal in the sense that rrvf 0 for every 

i G Nnt and j G Nn„.. 

Assumption 2 means that G is a block-diagonal matrix with the rank one 

matrices Mi from (3) along the diagonal 


Ew = 

/Ml 

0 

0 .. 
M2 .. 

• ° ^ 
0 

with w = 

/ mi \ 


U 

0 .. 

Mn,, j 


\WnJ 


and the other matrices in (lb) are given as 



fcA 


(dA 


/ ei \ 

Ex = 

1 

: Eu — 

: 

, e = 

; 


[Cnj 


\Dnb/ 


\^nb / 


The solution set of the complementarity problem (lb) is simply the cartesian product of the 
solution sets of (3). 

The requirement that mf be nonzero is without loss of generality and serves to avoid the 
presence of spurious degenerate complementarity variables. Assume = 0 and cl^’'^x -f 
Df’ u -t- ef = 0 in (3), then wf is unbounded above and completely independent from 
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all other complementarity constraints. This means its corresponding column in has to he 
identically zero to satisfy Assumption 1. Hence, could simply he eliminated from the 
model. 

Assumption 2 obviously limits what LC models (1) we can consider, e.g. a consequence of 
Assumption 2 is that = E^ A 0 which is generally not the case. We will prove in 
Section V that every continuous PWA system can he written as an LC model in the form (1) 
satisfying Assumption 2. Hence, while we have limited the number of LC models (1) we 
consider, we still cover all cases where equivalent continuous PWA dynamics exist. 

It is easy to prove the following result which provides an easy procedure to verify whether 
certain LC models satisfy Assumption 1: 

Lemma 1 (nullspace representation of Assumption 1). An LC model (1) with E^ = E^ A 0, 
e.g. one that satisfies Assumption 2, satisfies Assumption 1 if 

N{E^)CN{B^). 


Proof Use [9, Thm. 3.1.7]. □ 

Finally, we make an assumption on the existence of special solutions for the complementarity 
sub-problems (3): 

Assumption 3. For every fixed x G and u G M”" and every sub-problem (3) with i G 
there exists a j G such that 

+ ep'' < 0. 

This assumption ensures that for any (x, u) the solution to (3) is non-trivial, i.e. Wi / 0. 
Similarly to Assumption 2, Assumption 3 restricts the LC models (1) considered but will be 
satisfied by LC models obtained from an arbitrary continuous PWA model via our construction 
described in Section V. 

Failure of MFCQ for (2) means that the set of Lagrange multipliers is either unbounded 
(which can be a numerical problem) or may even be empty at a local optimum. This means 
it is generally not possible to characterize optimality for MPCCs using the classical KKT 
conditions. The main result of this paper relieves this problem for the particular optimal 
control MPCC (2): 

Theorem 1 (strong stationarity conditions). Let v* := (x*, u*, w*) be feasible for the optimal 
control problem (2) for an LC model satisfying Assumptions 1, 2, and 3. Then v* is locally 
optimal if and only if the classical KKT conditions in the sense of [5, Sec. 5.5.3], [2, Sec. 
5.1] for (2) admit a primal-dual solution pair. 

We will need a number of intermediate results to prove Theorem 1, so defer the proof to 
Section III. Theorem I is remarkable since it states that the classical KKT conditions are both 
necessary and sufficient local optimality conditions for the optimal control MPCC (2). This 
is in contrast to general MPCCs, for which the KKT conditions are sufficient for optimality. 
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but typically not necessary [50, Ex. 3]. We will present the KKT conditions from [5, Sec. 
5.5.3], [2, Sec. 5.1] for (2) in a concise manner in Section III. 


II. Structure of the Linear Complementarity Solutions 

For the proof of Theorem 1 we will need a number of results that characterize the solution set 
of the complementarity constraint (2c) under the assumptions made in the previous section. 

From Assumption 2 it is clear that we can focus on an individual complementarity sub¬ 
problem (3). Furthermore, for the structural results we derive below the actual values of the 
state X and control input u are of no consequence. Hence, to simplify the notation throughout 
this section we will consider the linear complementarity problem (LCP) 

0 < Mw + q _L m > 0, (4) 

where w G M G and q G M”. Note that we drop the sub-problem index i from (3) 

and assume x and u fixed to obtain 

q := CiX + DiU + e*. 

Moreover by Assumptions 2 and 3 

G N„: q^^ < 0 and M = mrrJ ^ 0 

for some m G 

Following the conventions of [9], we call the solution set of the LCP (4) for a fixed M and q 

SOL (g, M) := |t(; G w >0, Mw -|- (? > 0, w'^ {Mw + (?) = o| . 

From M = M~^ A 0 if follows fhaf SOL (g, M) is a convex polylope: 

Lemma 2 (LCP solulion sel). The solution set for the complementarity problem (4) is the 
convex polytope 

SOL (q, M) = {re G \ q^{w — w) = 0, {w — li)) = 0} 
where w is any solution of (4). 

Proof The fad lhal SOL ((?, M) is a convex polylope follows from MAO and [9, Thm. 
3.1.7]. We can simplify Ihe nolalion from [9] by noling lhal N{M) = N (m^). □ 


An immediate consequence of Lemma 2 is lhal Ihe left-hand term {Mw + q) in (4) has Ihe 
same value for every w G SOL {q, M) [9, Thm. 3.4.4]. 

For a fixed w G SOL {q, M) we inlroduce Ihe index sels 


a{w) 

:= {i G Nn 

Ml* 

(3{w) 

:= {i G Nn 

Ml* 

y{w) 

:= {i G Nn 

Ml* 


'^w -h (?t*l = 0, > 0}, 

'^w + (?t*l = 0, = 0}, 

-h (?t*l > 0, = 0}, 


( 5 ) 
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which partition N^. Complementarity constraints for which i G /3(m) are called biactive, and 
the set j3{w) is the biactive set. In case I3{w) is empty, w is called nondegenerate [9]. 

Lemma 3 (nondegenerate tangent vectors). Let w G SOL {q, M) solve the LCP (4) with 

/3(m) / 0, 

i.e. assume there is at least one biactive complementarity constraint. For every j G /?(re) 
there exists another w G SOL (q, M) such that 

a{w) = a{w) U {j} and /3{w) = I3{w) \ {j}. 


Proof. We will construct a vector 5 G such that w = {w + 5) satisfies the conditions of 
the lemma. From Lemma 2, any such 5 must satisfy re = (m + (5) > 0 and 5 G N {[m 
in addition to the conditions on a{w) and j3{w). 


Assumption 3 ensures that a{w) is always nonempty. Select an arbitrary i G a{w) and 
j G l3{w), and restrict 5 to he zero aside from the elements jW and We will choose the 
non-zero entries of 5 such that 


Sli 


< 


and > 0, 


( 6 ) 


which is sufficient to ensure both the non-negativity of w = w + 5 and the required conditions 

on a{w) and j3{w). 


Note that gW 
amounts to 


—mW(m^re) (analogously for so that the nullspace condition on 5 


—mW (mJw) 


SU] 


= 0 . 


(7) 


Recalling that Assumption 2 ensures that the vector m is element-wise non-zero, it is then 
sufficient to choose 





which satisfies both (6) and (7). 







□ 


Lemma 3 proves that in the presence of biactive complementarity constraints it is possible to 
find a tangential direction to SOL {q, M) such that exactly one index moves from /3 to a. This 
technical result is necessary to later prove our results on the tangent cone of the constraints 
in (2). Repeated application of this result provides the following corollary: 

Corollary 4 (strict complementarity). ^SOL {q, M) is non-empty, then there exists aw^ SOL {q, M) 
with (3{w) = 0. 



III. Optimality Conditions for the Optimal Control Problem 


In this section we will derive necessary and sufficient optimality conditions for the optimal 
control problem (2). The problem (2) is a specific case of the following general class of ajfine 
MPCCs: 


min J{v) (8a) 

V 

S.t. FinV + fin < 0, FeqV + feq = 0, (8b) 

Gv + g > 0, Hv + h>0, (8c) 

{Gv + g)'^ {Hv + h) = 0 (8d) 

Here, v G is the decision variable in (8), fin G feq G g,h £ M™, and all other 
quantities have compatible dimension. For the optimal control problem (2) with v := (u, x, w), 
we have 

G = (In ® Eu In ® Ex In ® Em) , g = In ® e, (9a) 

= (0 0 Inu^) , /i = 0. (9b) 


The other matrices can be constructed analogously. In the sequel we will require a num¬ 
ber of concepts from the MPCC literature, and will make reference to the general affine 
MPCC (8) to this end. The optimization problem (8) is nonconvex due to the complemen¬ 
tarity constraint (8d). Note that there are many slightly different but equivalent formulations 
of (8), in particular regarding the complementarity constraint. We will use the formulation (8) 
throughout the paper, which is without loss of generality (see [42]). 


A. Preliminaries: stationarity conditions for MPCCs 


We introduce a number of index sets relating to the active (complementarity) constraints for 
a given feasible v*\ 


{i G Nm 

o' 

II 

+ 

* 

tb 


> 0 } 

{i G Nm 

o' 

II 

+ 

* 

1 b 

H^E^y* + h\ 

* 1 = 0 } 

{i G Nm 

1 G^E]y* + g[i] > 0 , 

H^E^y* + h\ 

* 1 = 0 } 

{i G Np 1 

fW = o|. 

in 'Jin J 




( 10 ) 


Due to the analogy between these definitions and (5) the same letters are generally used in the 
literature to identify the singularly active and biactive complementarity constraints. Note that 
we omit the dependence of these sets on the feasible v* in the interest of a simpler notation 
whenever it is unambiguous. When /3 = 0 we say strict complementarity holds [42]. 


We define the Lagrangian function for (8) according to [2, 5] as 

C{v, r), /I, VG, Ih, 0 := J(v) -f fj~^ {EinV + fin) + {FeqV + feq) 

— og {Gv F g) — ijj {Hv F h) F ^ {Gv -|- g) {Hv F h ), 
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where we introduced KKT multiplier variables f/ G fi G uq G uh £ and 
^ G M for the constraints in (8). The classical KKT conditions for (8) are then given as 
follows [5, Sec. 5.5.3], [2, Sec. 5.1]: 


VvC{v, fj, fi, UG, 1) = 0, 


> 0 

for 

'1 £ 

fjlH 

= 0 

for 

t 0 ^irij 

= 0 

for 

i G 7, 

At] 

= 0 

for 

i G a, 

> 0 

for 

i G a U /), 

At] 

> 0 

for 

i G /3 U 7, 


( 11 ) 


A set of primal variables v* G feasible in (8) and KKT multiplier variables (also called 
Lagrange multipliers) that together satisfy the KKT conditions (11) are also called a primal- 
dual solution pair for (8). As discussed in the introduction the KKT conditions (11) may not 
be necessary optimality conditions for (8) due to failure of MFCQ [8]. A number of alternative 
stationarity conditions to the classical KKT conditions have therefore been introduced. We will 
limit the discussion to the concepts necessary for our control context; for a comprehensive 
discussion of weaker stationarity conditions that are applicable for more general MPCCs 
see [22, 55]. 


The strongest optimality conditions that can generally be expected to hold at an optimal 
point [37] are the so-called M(ordukhovich)-s,tntionwy conditions: 

Definition 1 (M-stationarity). A feasible point v* G M”” of (8) is called M-stationary if there 
exist rj G ML, p G vg £ vh £ such that 

VJ{v*) + fIp + Fjgp - G'^ug - H^vh = 0, 

>0 for i G Tin, r/W =0 for i 0 

fil [il (12) 

Uq = 0 for i G 7 , for i G a, 

{^'g ^ 0 ^ 0) V = 0 for i £ (3. 


The quantities ry, ug, and i/h are commonly referred to as Lagrange multipliers, although they 
are distinct from the classical Lagrange multipliers appearing in the KKT conditions (11). 
To avoid confusion we refer to the former here as the MPCC multipliers and use the term 
KKT multipliers for the latter. The constraints on the signs of these MPCC multipliers are 
different, i.e. the multipliers for the singularly active complementarity constraints, Vq for 

i G a and z/K for i G 7, are allowed to be negative. The same holds for one of the multipliers 

[ 2 ! [ 7 ] 

corresponding to the biactive complementarity constraints, Vq and for i G /?, as long 
as the other one is equal to zero. While it can be shown that the M-stationarity conditions 
correspond to non-smooth KKT conditions for an equivalent formulation of (8) [55], it is 
possible that an M-stationary local minimum v* of (8) does not admit a primal-dual solution 
to the classical KKT conditions (11) [50, Ex. 3]. 

The so-called strong stationarity conditions are stronger than those of M-stationarity and have 
a close relation to the KKT conditions of (8): 
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Definition 2 (S-stationarity). A feasible point v* G M”” of (8) is strongly (or S-)stationary if 
there exist p G p, G oq G M™, and vh G M™' such that 


VJiv*) + + Fj^p-G'aG-H'aH = 0, 

7?W >0 for i G Tin, p^"'^ =0 for i 0 Tin, 


=0 ^ 7 , 


'H 


= 0 for i G a, 


(13) 


^ > 0 A 1/K > 0 /or i £ [3. 


'G 


The only difference between M- and S-stationarity lies in the conditions on the biactive 
multipliers. Hence, the two stationarity conditions (and in fact most other MPCC stationarity 
conditions) collapse into S-stationarity in the case of strict complementarity (/? = 0). 

A feasible point v* of (8) is an S-stationary point if and only if there exist KKT multipliers 
satisfying the classical KKT conditions (11) for (8) at the same point [20, Prop. 4.2]. However, 
the classical KKT multipliers in (11) are distinct from the MPCC multipliers certifying 
S-stationarity in (13), e.g. the MPCC multipliers characterizing an S-stationary point do 
not include a multiplier for the orthogonality constraint (8d). It can be shown that MPCC 
multipliers vg and certifying strong stationarity of v* can be computed from KKT 
multipliers satisfying (11) as 

VG = >^G - i{Hv*+ h), VH = >^H - i{Gv* + g). (14) 

While Xg and Xh must be nonnegative because they originate from the classical KKT 
conditions, individual components of vg and vh might be negative as indicated in Definition 2. 
Derivation of KKT multipliers from given MPCC multipliers certifying S-stationarity is 
similarly possible [20]. 

While M-stationarity is a necessary optimality condition for the affine MPCC (8) [55], 
additional conditions on the MPCC multipliers are required to make M-stationarity a sufficient 
condition for local optimality, e.g. [55, Thm. 2.3]. S-stationarity, on the other hand, is a 
sufficient optimality condition for the affine MPCC (8), but is not a necessary optimality 
condition in general [50, Ex. 3]. We show, however, that in the particular problem (2) 
considered here the S-stationarity conditions are both necessary and sufficient for optimality 
as stated in Theorem 1. 


Our proof of Theorem 1 is based on a regularity property of the linearized tangent cone 
q-hnf^v*'^ to the Constraints of (8) at a feasible v* [23] defined as 


Vi G . 

Ft'd<0 


Feqd = 0 

yi G a: 

G^F]d = Q 

Vi G 7: 

H^F]d=Q 

Vi G /): 

G^F]d > 0 

Vi G /): 

> 0 
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and of the related MPEC-linearized tangent cone (n*) [20, 55]: 


Ti 


iin 
MPECK 


[v*) := n {d G M”” Vi G /3: = o} (15) 


The tangent cone T{v*) to the feasible set V of (8) is a closed cone [45, 49] and defined as 
Tiv*) := (dGM 


Vh — V 

3tk G M+, 'Ufc E V: lim tfc = 0 A lim - = d 

k^oo k^oo t]^ 


It holds that r(r;*) C Tj(tpEciv*) ^ [22] and, hence, (r'*^(r;*))° C {T^^pEciv*))^ C 

(T{v*))°, where C° is the polar cone of C [49, Sec. 6.E]. Also note that is a convex 

polyhedral cone while T{v*) and Tj}^Ec('*^*) generally nonconvex. 


B. Necessary and sufficient optimality conditions for the optimal control problem 

We can now prove an important geometric property of the feasible set of the optimal control 
problem (2). Specifically, that it satisfies the intersection property [21]: 

Lemma 5 (intersection property). The feasible set of (2) for an LC model satisfying As¬ 
sumptions 1, 2, and 3 satisfies the intersection property at every feasible v := (u, x, w), 

(A}?>i5c(v))° = (r'*"(v))°. (16) 

Proof If /3 = 0 then the result is obvious, so we assume j3 $ throughout and will use 
Lemma 1 from [45] to prove that equality holds in (16). To this end we will for every j G /? 
find a vector d G M”” such that 

Feqd = 0, An A = 0 Vi G Tin (17a) 

= 0 Vi G a U /?, = 0 Vi G 7 U /3 \ {j}, H^’^d > 0, (17b) 

where we used the more compact notation of (8). It can be seen that the d we construct must 
be a local recession vector to the feasible set of (2) such that the biactive complementarity 
constraint j becomes singularly active, i.e. for a sufficiently small e we have 

a{v + ed) = a{v) U {j} and /3{v + ed) = j3{v) \ {j}. 

The similarity with Lemma 3 is not coincidental, we will in fact use it in the proof. 

Let the recession vector d := (du, dx, d^) be such that du = 0 and dx = 0. Comparing (2) 
with (17) it is easy to see that (17) in that case decouples over the horizon. The complemen¬ 
tarity constraint j we are considering is part of (2c) for a particular timestep k and we can 
set all components of d^ := {dwo^ ■ ■ ■ ffiwN-i) riot corresponding to k equal to zero. 

We have now reduced the conditions (17) to finding a special recession direction for a 
single LCP (2c). Due to Assumption 2 this problem decouples even further into the linear 
complementarity sub-problems (3) in the variables *. Hence, by Assumption 3 and Lemma 3 
there exists a direction d satisfying (17b). 
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From the proof of Lemma 3 we can see that this G J\f{Ew) and, hence, hy Lemma 1 
also dwk £ Any additional linear inequality constraints present in the optimal control 

problem (2) would constrain only states and control inputs and the corresponding components 
of the recession vector d are assumed zero. 

The argument above proves that we can construct for any j £ f3 a recession vector satis¬ 
fying (17). This is equivalent to condition (Ah) in Lemma 1 of [45] with I3gHi = 0 and 
/?GH 2 = /3- By the same Lemma and Theorem 1 in the same reference it follows that the 
intersection property holds. □ 

With the intersection property established for (2) we are now in a position to prove Theorem 1: 

Proof of Theorem 1: Consider the tangent cone T (v*) to the constraints in (2) at the feasible 
V* from the theorem (for details see [21, 49]). We have, e.g. from [22], 

r(v*) c r^^pEci^*) c r'*"(v*). ds) 

Since (2) is an affine MPCC the MPEC-Abadie constraint qualification holds at every feasible 
V [23, Thm. 3.2] and in particular at v*, which means T(v*) = T^pEci'^*) [23, Def. 3.1]. 
Substituting this and considering the polar cones in (18) we obtain 

(T(V))° = (7j;?,£c(V))° 2 (G“(v*))°. 

From Lemma 5 we know that the intersection property holds for (2) and, hence, the inclusion 
is actually an equality and we finally have 

This is fhe so-called Guignard constraint qualification. If follows from [21, Thm. 16] fhat 
S-sfafionarify, in fhe sense of Definilion 2, is a necessary opfimalify condition for (2). We 
also know from [55, Thm. 2.3] fhaf S-sfafionarity is a sufficienl opfimalify condition for (2). 
S-sfationarify is in facf equivalenf fo fhe classical KKT-condifions for (2) [20, Prop. 4.2]. This 
complefes fhe proof. □ 

Nofe fhaf fhe conclusion of Theorem 1 also holds wifh addifional linear constraints on the 
states and inputs, because the recession direction d constructed in the proof of Lemma 5 has 
no components in the x and u directions, i.e. the intersection property will still hold. 


C. Stronger conditions for global and isolated optima 

We now present a number of stronger sufficient conditions that can be used to verify whether 
a given locally optimal solution v* ;= (u*, x*, w*) is, in fact, globally optimal or an isolated 
minimizer. The latter means that there is no other locally optimal point within a neighborhood 
around v*. To present sufficient conditions for an S-stationary solution v* of (2) to be globally 
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optimal, we introduce MPCC multipliers Hk G for (2b) and for (2c) to 

present the S-stationarity conditions for v*: 

— {x*n) + /tat-i = 0 (19a) 

dih 

\/k = l,...,N -1: ^{xl,ul) + Hk-i-A^f^k-Eji^k = 0 (19b) 

dll. 

yk = 0,...,N-l: ^(^xl,ul)-Bjfik-Ejuk = 0 (19c) 

yk = 0,...,N -1: - Afc = 0 (19d) 

We omit the (inconsequential) conditions related to the derivatives with respect to xq. Since 
the complementarity constraints in ( 2 c) naturally decouple over the prediction horizon for 
the given S-stationary v*, we can consider index sets fik, and 7 ^ analogously defined as 
in (10) for each stagewise LCP. With that, we have the following conditions for the MPCC 
multipliers corresponding to the complementarity constraints ( 2 c): 

Vi G afc: = 0 , Vi G 7 ^: 4*^ = 0 , (19e) 

ViG/3fc: (19f) 

We can now state a sufficient condition for an S-stationary solution of (2) to be globally 
optimal. 

Theorem 2 (sufficient global optimality condition). Let v* = (u*,x*, w*) be an S-stationary 
point for (2), i.e. there exist MPCC multipliers pk ^ Uk G M”™, and Xk G such 
that (19) holds. If additionally for all k G Ntv-i 

Ok > 0 and Xk > 0, (20) 

then V* is a globally optimal solution to (2). 

Proof Since all constraint functions in (2) are affine and we assumed convex stage- and 
terminal cost functions Ik the result follows from [55, Thm. 2.3]. □ 


In addition to the first-order optimality conditions from Theorems 1 and 2 we can use second- 
order sufficient conditions to identify an isolated minimizer of (2). These pose conditions on 
the critical directions at a stationary point v* of ( 8 ) to ensure that no feasible non-ascent 
direction exists. To this end, we introduce the so-called MPEC critical cone based on the 
MPEC-linearized tangent cone from (15): 


Cmpec{v*) :— Vi G M"" 


VJ(u*)^(i< 0} 


With this we introduce the following second-order sufficient condition: 


Definition 3 (strong second-order sufficient condition). Given an M-stationary point v* G 
for ( 8 ) we say that the M-multiplier strong second-order sufficient condition (M-SSOSC) holds 
at V* if and only if 


yd(^CMPEc{v*)\{^}. d~^V^J{v*)d> 0. 
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Note that for (8), M-SSOSC does not depend on the values of any MPCC multipliers hut 
instead only on the critical directions. 

For a comprehensive discussion of this and other second-order conditions for general MPCCs 
we refer to [27]. Note that we have simplified substantially the definitions from [27] since we 
are only considering the affine MPCC (8). We can use the M-SSOSC from Definition 3 to 
identify an isolated minimizer v* of (2), i.e. the only local minimizer within a neighborhood 
around v*. 

Theorem 3 (isolated minimizer with unique complementarity variables). Let v* = (u*, x*, w*) 
be an S-stationary point for (2) and let M-SSOSC hold at v*. Then v* is an isolated local 
minimizer, i.e. there exists an e > 0 such that 

N-l N-1 

||v-v*||<e hixl,ul) < iNjxN) + ik{xk,Uk). 

k=0 k=0 

Furthermore, w* are the only complementarity variables solving the TCP (2c) for the fixed 
u* and X*. 

Proof. The fact that v* is an isolated minimizer follows from [27, Cor. 4.3]. 

Assume there exists another w such that v := (u*,x*,w) is feasible for (2). Since the 
LCP (2c) decouples over the horizon by Lemma 2 the set of feasible w is a convex polytope, 
hence w can be arbitrarily close to w*. Since w does not enter the cost-function (2a) we 
have J(v) = J(v*) and reached a contradiction. □ 

IV. Optimality Conditions for Optimal Control Inputs 

The proof of Theorem 3 hints at a possible challenge when solving the optimal control 
problem (2) that, to the authors’ knowledge, has not yet been considered in the literature: 
we are interested in finding optimal control input trajectories u* and the corresponding state 
trajectory x*, but the optimal control problem also includes the optimal complementarity 
variables w*. Fundamentally, the problem we want to solve is the hybrid optimal control 
problem 

7V-1 

min ^Ar(xAr)-h 4(xfc,Ufc) (21a) 

U,X ^ 

s.t. V/c = 0,..., - 1: Xk+i = f{xk,Uk), (21b) 

where /: x M”" —)• describes the LC dynamics (1). The problem (2) is an in¬ 

stance of (21) that explicitly considers complementarity variables w to facilitate theoretical 
derivations and numerical computations. 

It is conceivable that a local optimum v* = (u*, x*, w*) of (2), which is a member of the 
feasible set 

Q := {v = (u,x, w) I V satisfies the constraints in (2)} , 
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has no corresponding optimal solution (u*,x*) to ( 21 ) in the set 

U := {(u, x) I 3w: (u, x, w) G Q} , 

which is the projection of Q onto the state-input space and the feasible set of ( 21 ). 

We will next present an example to illustrate this point. To this end we introduce the set 

Ad(u, x) := {w I (u,x, w) G Q} 

of complementarity variables that are consistent in ( 2 ) with a given u and x. 

A. Example: non-optimal input trajectories can be S-stationary 
Consider the LC model 


'+ = + u;[2l _ 2, 



( 22 a) 

0 < re'^' -|- + x + u 

± 

reW > 0, 

( 22 b) 

0 < rr;[^] -h - 1 

± 

rr;[2] > 0, 

( 22 c) 


which is equivalent to the continuous PWA model 

x"*" = max{—(x + u + 2), —1}. 

It is easy to verify that (22) satisfies Assumptions 1-3, the details are omitted here for brevity. 

Consider the optimal control problem (2) with initial state xq = 0, prediction horizon = 1, 
stage cost £q{xo,uo) = and terminal cost ('i(xi) = ^xf. For a control input uq = —1 
we obtain xi = — 1 and an infinite number of admissible complementarity variables with 

M(-l, -1) = {wo G I + reP = l} . (23) 

Solving the S-stationarity conditions (13) for any choice of wq G Af(—1,—1) yields the 
unique MPCC multiplier candidates 

Eo = 1, 1^0 = , Ao = . (24) 

To analyze optimality we have to distinguish three cases: 

Case 1: = 1, =0 a = {!}, /3 = { 2 }, 7 = 0 

Case 2: > 0, >0 ^ a = {1, 2 }, /3 = 0, 7 = 0 

Case 3: rug^^ = 0, = 1 => a = {2}, /3 = {!}, 7 = 0 

Case 2 is the easiest to analyze since the biactive set jj is empty. Hence, for any ref' > 0 

and ref' > 0 , the point v = (—1,—l,reo) is S-stationary and therefore a local optimum. 

\ 2 ] [21 

For case 1, we see that the MPCC multipliers Uq and Ag corresponding to the biactive 
complementarity constraint 2 are both zero, hence, in this case we are also at an S-stationary 
local optimum. Case 3, on the other hand, does not satisfy the conditions for S-stationarity (or 
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alternative sufficient optimality conditions, e.g. from [27, 55]) since < 0. From Theorem 1 
it follows that v* = (—1, —1,0,1) is not locally optimal since it is only M-stationary and not 
S-stationary. 

Figure 1 shows the cost function values J(x, u) in (2a) for varying values of uq (which is 
also the value function of ( 21 )) together with the corresponding admissible values of 
Corresponding values for Wq and xi follow from the dynamics (22). The three cases discussed 
above are marked in green along the dashed red line indicating the feasible set of the optimal 
control problem. 



-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 

Uq 


Fig. 1. Objective function value of optimal control problem (2) for example system (22) and varying control 
inputs together with admissible values for 

From the value function J in Figure 1 it is clear that a control input of uq = — 1 should not 
be considered locally optimal from a control perspective and will not correspond to a local 
optimum of ( 21 ) in (7 because ft = — 1 -F e for an arbitrarily small e > 0 is feasible and has a 
lower objective function value. At the same time, it illustrates how the local minima for (2) 
in cases 1 and 2 come about: in a small neighborhood around either of these points, a control 
input rio > — 1 is not feasible because this would require a significant change in . In other 
words, while 

||(u*,x*) - (u,x)|| 

might be arbitrarily small, 

||(u*,x*, w*) - (u,x,w)|| > e Vw e Ad(u,x), 
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for some finite e > 0. Hence, there exist values for wo such that v* = {uq, xi, wq) is locally 
optimal for (2) with uq = —1. For = 0 this suddenly changes because there exist 
tangential descent directions pointing inside the feasible set. At the same time, for any such 
V* G Q the corresponding point (no,xi) G is not a local optimum for ( 21 ). 


B. Necessary and sufficient conditions for optimal control trajectories 

The preceding example highlights the necessity of adapting Theorem 1 to the control setting. 
For control applications, one is primarily interested in the values of u* and x* and can 
disregard the actual value of w* (provided it is primal feasible). To that end, we introduce 
the notion of a locally optimal control trajectory. 

Definition 4 (locally optimal control inputs). A control input trajectory u* := 
is called locally optimal with respect to (2) if there exists an e > 0 such that 

N-l N-l 

^n{x*m) + ^ ik{x*k,ul) < In{xn) + ^ ik{Xk,Uk) 

k=0 k=0 

for all control trajectories u ;= {uo,ui,... ,un-i) with ||u — u*|| < e. Here, x* and x are 
the state trajectories resulting from applying control inputs u* and u, respectively. 

Obviously, a locally optimal control trajectory corresponds to local optima of (2) and (21). The 
example in Section IV-A illustrates that the reverse does not hold: there may exist local optima 
for (2) that do not correspond to locally optimal control trajectories in (21). We immediately 
obtain necessary and sufficient conditions for a state and input trajectory (x*, u*) to be locally 
optimal. 

Theorem 4 (optimality conditions for input trajectories). Given an input trajectory u* G 
corresponding state trajectory x* G for an LC model (1) satisfying 

Assumptions 1—3, u* is locally optimal per Definition 4 if and only if v* := (u*,x*,w) is 
S-stationary for all w G Af(u*,x*). 

Proof. The input trajectory u* being locally optimal is equivalent to v* being a local optimum 
of (2) for all w G Af(u*,x*). The result follows from Theorem 1. □ 

A number of sufficient conditions for an input trajectory u* to be locally optimal can be 
derived which are very easy to check. The proofs of these results are straightforward. The 
first Corollary mirrors the sufficient conditions in Theorem 2 for a point to be globally optimal 
in ( 2 ). 

Corollary 6 (globally optimal control trajectories). Let v* = (u*,x*,w*) be an S-stationary 
point for (2) and let the MPCC multipliers satisfy the conditions (20) in Theorem 2. Then u* 
is a globally optimal input trajectory. 


Specializing Theorem 3 to optimal control inputs yields the next result. 
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Corollary 7 (isolated control trajectories are optimal). Let v* = (u*,x*,w*) be an S- 
stationary point for (2) and let M-SSOSC hold at v*. Then u* is an isolated locally optimal 
input trajectory. 

While isolated local minimizers always have unique complementarity variables (see The¬ 
orem 3) other (non-isolated) local minimizers can have the same property. Hence, from 
Definition 4 follows Corollary 8. 

Corollary 8 (unique complementarity variables imply optimality). Let v* = (u*,x*,w*) 
be an S-stationary solution for (2) and let Af(u*,x*) be a singleton, i.e. w* is the unique 
complementarity variable trajectory consistent with u* and x*. Then u* is a locally optimal 
input trajectory. 

Theorem 3.1.7(b) in [9] can be used to determine whether the given complementarity variable 
trajectory w* is unique using the determinant of a particular submatrix of Ln (g) Ew. From 
Corollary 4 we know that /3 = 0 is a necessary condition for Ad(u*,x*) to be a singleton. 

V. Inverse Optimization Modeling of Piecewise Affine Systems 

In this section we will show how all the results from the previous sections can be applied to 
hybrid dynamical systems in continuous piecewise affine (PWA) form. Their state dynamics 
are given as 

x'^ = f{x, u) = AiX + BiU + Ci for {x, u) G Oj, (25) 

where x G is the system state and u G the control input. The Ur regions D, C 
xM”“ form a partition of the domain Ll of the PWA system (25), i.e. int(Oj)nint(Dj) = 0 
for i / j and [fflli O* = 0. We assume the system dynamics to be continuous across region 
boundaries, i.e. 

AiX + BiU + Ci = AjX + BjU + Cj 

for all X G Dj n Dj. 

The finite horizon optimal control problem for PWA system (25) is given as 

N-l 

min iNixw) + h{xk,Uk) (26a) 

U,X ^ 

^=0 

s.t. yk = 0,... ,N - 1: Xk+i = AiXk + BiUk + Q for (x^, Uk) G D*, (26b) 

Here, u := (uq, • •., ua-i) and x := (xq, • • •, xa) are the input and state trajectory from an 
initial state xq as in (2). These problems are typically solved with mixed-integer program¬ 
ming (MIP) approaches based on an MLD reformulation of the PWA dynamics [1]. Other 
solution approaches based on nonlinear programming [14], the solution of a series of linear 
programs [15], or the alternating direction method of multipliers [24] have also been proposed. 

While the equivalence between PWA models (25) and LC models (1) has been known in the 
literature for quite some time, the LC models resulting from the derivations in [30] will not 
satisfy Assumptions 2 and 3 and require an MLD model as an intermediate step. The authors 
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in [30] point out that there will generally he a multitude of LC models (1) corresponding to a 
given PWA model (25). Which of these LC models is “best” will depend on the application. 

Recent results from inverse optimization [31, 35] provide a direct link between PWA mod¬ 
els (25) and LC models (1). This approach is based on the difference of convex functions [38] 
and in most cases of interest provides very compact LC models [31, Lem. 4]. We can show 
that following this approach will lead to an LC model that satisfies Assumptions 1-3. 

We will later make use of the results in [35] to represent the PWA system (25) as an optimizing 
process. To this end we will require the following result on the representation of a continuous 
PWA function —)■ M” with only convex component functions —)■ M as the 

optimal solution to a parametric quadratic program: 

Lemma 9 (convex PWA function as solution to PQP). Let V’: — )■ M” be continuous PWA 

and such that —)■ M is convex for all i G N^. Then 

1 11 — 112 

V’(p) G arg min - \\y — 'f{p)\\n sh. y > 'f{p) (27a) 

j/sR* 2 ^ 

for any diagonal matrix Q y 0 and an affine function —)• with 

f{p) < 'il){p) \/p. (27b) 

Furthermore, the arg min is a singleton, i.e. fip) is the unique optimizer for (27a). 

Proof Since Q 0 is diagonal (27a) decouples into scalar optimization problems in yW as 

follows: 2 

G arg min fyW — s.t. 

j/HeR 2 V / 

It is easy to see that this is simply the projection of VjW(p) onto the (convex) epigraph of 
The result follows immediately from (27b). □ 

Lemma 9 is a variant of [33, Lem. 2] which considered a parametric linear program in (27a). 
Note that we can always find an affine function satisfying (27b) due to the convex component 
functions of ijj, cf. [49, Thm. 8.13]. 

We will adapt the procedure outlined in [35] to derive a so-called inverse optimization model 
for a given PWA system (25). To this end we use the fact that scalar-valued continuous PWA 
functions can be written as the difference of two convex PWA functions [33, 40]. For the 
PWA dynamics (25) this means f{x,u) = ip{x,u) — (p{x,u) where 

ip{x, u) = Ayjx + Byju + Cyj for (x, u) £ Llj, (28a) 

4>{x, u) = A^^kX + B^^kU + for {x, u) G (28b) 

are continuous PWA functions over the domain of the PWA system and the component 

functions V'W, x —)■ M are convex for all i G The reader is referred to [31, 

33] for details on the computation of the convex decomposition including the matrices Ayj, 
Byj et alia. Note that this convex decomposition of / is not unique. The regions and 
Lip also form partitions of the domain for j G and k G Nnj. They may differ from 
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the original regions flj, unless the original regions form a so-called regular partition [35, 
Lem. 2]. The reader is referred to [11] for a comprehensive treatment of regular partitions; 
typical examples are Delaunay triangulations, Voronoi diagrams, and related partitions. 

We can apply Lemma 9 to the PWA functions ip and (p in (28) that make up the PWA 
dynamics (25) to obtain the following inverse optimization model of (25): 


x~^ = y — z (29a) 

y G arg min ^ lly - ^(x, u) 11^ (29h) 

s.t. y > Ay^jX + Byju + Cyj Mj G (29c) 

2 G arg min -\\z — p{x,u)\^„ (29d) 

s.t. z > A^^kX + B^^kU + V/c G Nn- (29e) 


From Lemma 9 it follows that Qy,Qz G are arbitrary positive definite diagonal 

matrices. The affine functions in (29b) and p in (29d) are given as 

ip{x, u) = A^x + B^u + c^, 
p{x, u) = A^x + B^u + c^, 

and musf safisfy (27b) wifh respecf fo V’ and cp, respectively. 

If follows direcfly from fhe consfrucfion fhaf fhe inverse optimization model uses n = 2nx 
decision variables and a sfrictly convex cosf function. For alternafive, possibly more compacf, 
construcfions of an inverse opfimization model of (25) see [34, 35, 44]. An immediately 
obvious simplification can be performed in case a component function /W of the original 
PWA dynamics is already convex or concave. 

Recall that (29) represents a remodeling of the PWA dynamics in (26b) in terms of a parametric 
optimization problem. To facilitate the inclusion of this model into our overall optimal control 
problem (26), we can rewrite (29) in terms of its KKT conditions, yielding the following linear 


complementarity model for our system dynamics: 

x^ = y — z, (30a) 

nl 

0 = Qy{y - {A^x + B^u +c^)) (30b) 

i=l 

0 = Qz{z - {A^x + B^u +c^)) (30c) 

j=i 

Vf G : 0 ^ y — (^Ay^ix By^iU -\~ Cy^i) -L Xi A 0, (30d) 

Vj G : 0 < z — {Azjx + Bzju + c^j) _L Oj > 0. (30e) 


To avoid confusion with the KKT multipliers and MPCC multipliers introduced in Section III, 
we will call the A* and 9j in (30) internal multipliers, i.e. multipliers for the lower-level part 
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of the optimal control problem (26) when the dynamics have been reformulated as the KKT 
conditions of a parametric QP. They correspond to the complementarity variable w in (1). 

Note that the reformulation (30) of the PWA system (25) is an LC model (1) with a generalized 
cone complementarity. It can be shown that all results in Section III also hold for the 
formulation (30) with the additional auxiliary variables y and z [32]. Alternatively, one could 
eliminate these auxiliary variables to obtain an equivalent, more compact complementarity 


representation that exactly matches (1): 

= {A^ - A^) X + - B^) u + c^- c^ + Q~^ Qz^ ^ Gj (31a) 

i=l j=l 

wx^ — Qy ^ ^ Aj A (^A^ Ay i'j X A {^Bqp By i'j u A (31b) 

i=i 

0 < w\. _L Ai > 0 Vi G (31c) 

A 

WQj = Qz^ 'y y 3" (^<i> “ ® V ~ ^z,j) u A C(f, — Czj (3 Id) 

i=l 

0 < W0^ A % > 0 Vj G (31e) 


What is left to do is to show that it also satisfies the assumptions made at the start of the 
paper. To that end, we make the following non-restrictive assumption on the functions 4) and 
4>: 

^{x,u) < max AyjxAByjuACyj V(x,u) G (32a) 

’ ’ ’ 

d>(x, u) < max Az iX A Bz iU A Cz i V(x, u) G 12 (32b) 

ieN„^ 

While (32) is stronger than (27b), it can always be satisfied for a given PWA system, e.g. by 
choosing any index j G and i G and setting 

r/!(x, u) = Ay^jX A ByjU A Cyj — y, (33a) 

(j){x, u) = A^^iX A B:z,iU A Cz,i - C, (33b) 

with arbitrary y,C, >0. With this we can prove that our results from Section III also apply 
to properly remodeled PWA systems: 

Theorem 5. Any LC model in the form (31) satisfying inequalities (32) also satisfies Assump¬ 
tions 1, 2, and 3. 

Proof For a fixed x G and u G M”" the conditions (30b) and (30d) form the KKT system 
of the strictly convex PQP (29b)-(29c) which has a unique minimizer, i.e. y is independent of 
the values of the internal multipliers A^. An analogous argument holds for the 9j in (31d)-(31e) 
and 2 in (29d)-(29e), hence (31) satisfies Assumption 1. 

Careful examination of (31) reveals that the internal multipliers always appear as their sum 
over the regions of the convex decomposition. Accordingly, coupling between the internal 
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multipliers occurs only for corresponding components. For any i G we can collect all 
internal multipliers influencing that particular component of x'^ in the following LCP: 


1 




V/c G : 0 < ^ 

Vj/ j=i 

Collecting terms appropriately results in a complementarity suh-prohlem (3) with 


Mi := 


Q 


[i,i\ 


ln“xn“ 


Q 






Ci := 


\ y,n^/ 


1 


1 1 


II 


II 


1 

W - 4 i) 

1 1 



An analogous argument holds for the LCP in 6, hence (31) satisfies Assumption 2. 

Finally, the component-wise satisfaction of (32) guarantees that for every i G there exists 
a /c G such that 

' 4 ^' - -+ { 4 "' - A:*) “+- 4 < »■ 

Again, we can show the same for the LCP in 0 and have proven that (31) satisfies Assump¬ 
tion 3. □ 


Theorem 5 finally proves our claim from the beginning of the paper — every continuous 
PWA model (25) can he rewritten as an equivalent LC model (1) satisfying Assumptions 1, 2, 
and 3. Accordingly, such LC models are very general and cover a wide range of applications. 
As will he seen in the next section, the slightly less compact hut sparser formulation (30) can 
he advantageous for computational purposes. 


The example in Section IV-A illustrates that the possibly infinite number of admissible internal 
multipliers for a given (x, u) can lead to spurious local minima in (2) that do not correspond 
to locally optimal control trajectories. It is easy to show that when writing a PWA system (25) 
as the LC model (31), different consistent internal multipliers only exist for a given (x, u) G 
when there exists an i G such that 


+ ^y,j^ + S'i = + ^y,k^ + 


for j 4 k in (28) (or analogously for z). This is naturally the case when (x, u) is on the 
boundary of two neighboring regions of the convex or concave part of the PWA dynamics. 
Accordingly, any control trajectory u that lands on such a region boundary will correspond 
to a local minimum in (2) with Af (u, x) not a singleton. It is conceivable that many of these 
local minima do not correspond to locally optimal control trajectories per Definition 4. While 
the simulations below indicate that this problem is not severe it is an issue that warrants future 
research. 
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VI. Numerical Results 

We first present a small toy example adapted from [33] that illustrates how PWA dynamics 
can he systematically transformed into an LC model. By construction, the resulting model 
will satisfy the conditions of Theorem 5 and, hence. Theorem 1. We then present extensive 
computational results on randomly generated PWA systems that show the possible henefits 
from reformulating PWA dynamics as descrihed in Section V. Using standard NLP solvers 
to solve the resulting MPCC (2) instead of dealing with a MIP formulation of (26) can he 
significantly faster without sacrificing much in terms of solution quality. 


A. Illustrative Example 

Consider the PWA dynamics / shown in red in the middle of Figure 2 which represents the 
dynamics of a hypothetical PWA system (25) with one state x and one control input u. Using 
the results from [33, 35], it can he decomposed into the two convex PWA functions i/; and cf) 
shown in the figure such that f{x,u) = 'ip{x,u) — (j){x,u). 


■ u) ■ f{x, u) u) 



Fig. 2. Example of PWA system dynamics / (in the middle) and its decomposition into a convex part ip (above) 
and a concave part —(p (below). 
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The explicit expressions for the functions V’: [—5, 5] x [—5, 5] —)■ M and (/>: [—5, 5] x [—5, 5] —)• 
M are given as follows: 

u) = max{3, x + 2, —x + 2, u + 2, —u + 2} 

(/>(x, u) = max{ 6 , x + tt + 2, x —tt + 2, 

—X — u + 2, —X + M + 2} 

Their difference yields the given dynamics /(•) = ip{■) — (/){■) which are omitted due to space 
limitations. 

Let us choose V'(x, u) := 1 and 0(x, tt) := x + u, which satisfy (27h) with respect to ip and 
(j), respectively. We can now use Lemma 9 to construct the following optimization problem 
in the two scalar decision variables zi and Z 2 from the expressions for ip and (p\ 

min ^ {zi - if + ]- {z 2 - (x + u)f (34) 
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(x, u) G [—5, 5] X [—5, 5] 

Due to its separability we can obtain the explicit solution to this optimization problems as 
zl{x,u) = ip{x,u) and z^ix^u) = (p{x,u). By construction, the PWA dynamics / shown in 
Figure 2 are recovered as /(x, u) = z^ (x, u) — (x, u) = ip{x, u) — (p{x, u). In other words, 

we have obtained an inverse optimization model (29) for the PWA dynamics /. 

It is straightforward to derive the LC model representations (30) and (31) from (34) by 
following the derivations in Section V. Because of our choice for ip and (p, the resulting LC 
model will by construction satisfy the conditions of Theorem 5. Hence, instead of solving 
the optimal control problem (26) for the original PWA dynamics we can instead solve the 
MPCC (2) and be assured that Theorem 1 holds. Accordingly, we have a good chance to 
solve the MPCC problem with standard NLP solvers which try to find solutions to the KKT 
conditions. 


B. Computational Experiments 

To investigate the possible computational benefits gained over traditional approaches from 
using an LC model (1) and solving the optimal control MPCC (2), we randomly generated 10 
different PWA models (25) with = 3 and = 1. The regions of these PWA models 
form a Delaunay-triangulation of the (bounded) system domain with 17 to 19 regions per 
system. We used the Multi-Parametric Toolbox 3.0 [36] and ECOS [18] for all geometric 
computations necessary to generate the systems. 
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To solve (26) we used the MLD reformulation from [1] and modeled the resulting MIP formu¬ 
lation with YALMIP [41]. For each prediction horizon A” G {2,4,, 10} we generated 50 
different initial states xq inside the domain of the system such that (26) was feasible, and 
solved the mixed-integer optimal control problem with Gurobi [28]. As the cost function we 
chose ^ ^ 

4(xfc,ttfc) = - (^\\xk\\l + \\uk\\fj and iN{xN) = ^ W^nWI ■ 


Each of the generated PWA systems was transformed into both the general LC form (30) and 
the compact form (31). The design parameters were chosen as Qy = Qz = I and and if) 
by shifting parts of the convex decomposition downwards as shown in (33). The resulting 
formulations of (2) were solved from the same initial states xq using the general purpose 
NLP solver IPOPT [54]. No special treatment of the complementarity constraints (2c) such 
as regularization [51] or penalization [47] was performed; we simply implemented them as 
scalar bilinear inequalities. 
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Fig. 3. Performance profile for A = 8. 

To get an impression of the performance that can be expected, Figure 3 shows a performance 
profile [17] of the computation times required to solve the 500 problem instances for N = 8. 
Define fhe performance rafio 
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where tp^s is the time it takes solver s to solve problem instance p. The performance profile 
then plots the function ps'. M —)■ [0,1] defined as 

Ps{t) := — \{p I Vp^s < t }\, 

Up 

which (for large numbers Up of problems p) is the probability that a performance ratio ^ 
is within a factor of r G M of the best possible ratio. 

It can be seen from Figure 3 that in over 95 % of all problem instances for = 8 the 
sparse LC formulation (30) was solved fastest. Additionally, for every problem instance the 
MIP approach is (in terms of computation time) beaten by either the sparse LC (30) or 
the compact LC (31) approach. It is worth noting that such a consistent performance of a 
general purpose NLP algorithm such as IPOPT is a direct result of the special formulation 
of the optimal control problem (2) and Theorem 1. It should not be expected when solving 
a general MPCC (8) because the KKT conditions (11) might not be satisfied at an optimum. 

In addition to the often slower solution times for the compact formulation (31), IPOPT 
encountered computational issues for 9 out of 2500 problem instances, e.g. reaching the 
maximum number of iterations, numerical issues etc. The sparse formulation (30) had no 
such issues and solved all instances to local optimality. While the sparse formulation (30) 
requires more decision variables in the optimal control problem (2) it also results in sparser 
constraint matrices whose structure may be more favorable to IPOPT. The issues of the 
compact formulation might be relieved by proper scaling or choice of the design parameters 
Qy, Qz, i’, and 0 in (29), but we will for the rest of this section only consider the sparse 
formulation (30). 

Figure 4 shows the best-case, worst-case, and median computation times necessary to solve 
the 500 problem instances for each prediction horizon for the sparse, general LC model (30) 
and the PWA model (25). The LC formulation (30) shows both much shorter computation 
times as well as gentler scaling for longer prediction horizons, e.g. the worst-case computation 
time for = 10 is two orders of magnitude larger for the MIP formulation. 

To evaluate solution quality we can consider the relative difference between the achieved 
objective function values 

J*NLP - J*MIP ^ 100 
’^MIP 

For 99.2 % of the 2500 problem instances over all prediction horizons this quantity is lower 
than 10 % and in 71.9 % of the cases the inverse optimization approach even yields a globally 
optimal solution. 

These results are affected by the choice of design parameters when constructing the inverse 
optimization model (29). For the results shown here we chose ^ and (p as shown in (33) with 
T] = ^ = 0.5 • 1, while with p = ^ = 10 ■ 1 less than 50 % of the problem instances are 
solved to global optimality. Referring back to our example in Section IV-A, the value of p 
corresponds to the (invariant) value of the sum of the internal multipliers in (23). Hence, a 
larger value of p corresponds to a larger set (in volume) of equivalent internal multipliers, 
i.e. cases 1 and 3 in Figure 1 move further apart. 
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Fig. 4. Computation times for different prediction horizons N. The center line is the median computation time 
over 500 problem instances for each prediction horizon and the error bars indicate the best and worst computation 
times, respectively. 


For 187 of the 249 problem instances whose optimal values are further than 1 % away 

from the globally optimal IPOPT returns control input trajectories unlp which for at 

least one time step land exactly on the boundary between two regions of the PWA system. 
The other 62 problem instances correspond to locally optimal input trajectories in the sense 
of Definition 4 that traverse the interior of the PWA system’s regions. To verify whether the 
boundary cases correspond to locally optimal trajectories we could use Theorem 4 (which 
requires a vertex enumeration for the set M.{\jlp!lP:'^nlp)) or one of its Corollaries. Due 
to the relatively short computation times for the optimal control MPCC it may be a viable 
strategy to solve (2) from different initial points and use the best result. This would reduce 
the chance of obtaining a solution on a region boundary and could even be implemented in 
parallel if the hardware allowed this. 

Using the LC formulation (30) instead of an MIP approach to the PWA model (25) presents 
a trade-off where a significantly shorter and more consistent solution time is bought with a 
potential decrease in solution quality. Simulation evidence indicates that this degradation is 
negligible and its severity within reasonable bounds in most cases, in particular because the 
optimal control problem (2) is typically solved in a receding horizon setting where only the 
first step of the control input trajectory u* would be applied before re-solving the problem at 
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the next time step [48]. 


VII. Conclusion 

In this paper we consider constrained optimal control problems for hybrid dynamical systems 
of linear complementarity type satisfying a number of structural assumptions. The optimal 
control problems are then mathematical programs with complementarity constraints which can 
be modeled as continuous nonlinear programs. Under the assumptions made, we can prove 
that the classical Karush-Kuhn-Tucker conditions are necessary and sufficient for optimality 
which is rarely the case for general MPCCs. 

Additionally, it is shown how continuous piecewise-affine systems can always be written as 
LC models satisfying our initial assumptions. This enables the treatment of control problems 
for a large and important class of hybrid systems as continuous NLPs using standard solution 
software. Numerical simulations illustrate the efficacy of this approach where the NLP can be 
solved in significantly shorter and more consistent time than a more traditional mixed-integer 
approach. The downside is a possible degradation in solution quality although this is often 
negligible since the NLP approach finds the global optimum in many cases. 

It should be investigated how much can be gained from using a solution method tailored 
to MPCCs. A number of relaxation and penalization methods have been suggested in the 
literature, cf. [37] for a recent overview, but they have to accommodate the fact that M- 
stationarity is the strongest optimality condition available for most MPCCs. A straightforward 
nonlinear programming approach to (2) can prove successful due to the strong stationarity 
results in this paper. Alternatively, an MPCC-method that guarantees convergence to S- 
stationary solutions could be used, cf. [39]. 
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